home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_400 / 420_02 / plot3d.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-01  |  38.4 KB  |  938 lines
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <math.h>
  4. #include <stdlib.h>
  5. #include <iostream.h>
  6. #include "titillat.h"
  7. #include "plot3d.h"
  8.  
  9. #ifndef TRUE
  10. #define TRUE -1
  11. #endif
  12. #ifndef FALSE
  13. #define FALSE 0
  14. #endif
  15.  
  16. typedef struct
  17.           {
  18.             int x;
  19.             int y;
  20.           } box_rec;
  21.  
  22. plot3d::plot3d()
  23.   {
  24.     prime_array_allocated=FALSE;
  25.     plot_prepared=FALSE;
  26.   }
  27.  
  28. plot3d::~plot3d()
  29.   {
  30.     if (prime_array_allocated)
  31.       delete prime_array;
  32.   }
  33.  
  34. int plot3d::prepare_plot(
  35.   double (*f)(double,double),
  36.   double x_min,
  37.   double x_max,
  38.   double y_min,
  39.   double y_max,
  40.   int    (*external_to_plot)(double,double),
  41.   int    (*red)(double,double),
  42.   int    x_division_count,
  43.   int    y_division_count,
  44.   double rotation_in_degrees,
  45.   double tilt_in_degrees,
  46.   double light_x,
  47.   double light_y,
  48.   double light_z)
  49. //      This function prepares a plot for generation.  If returns TRUE if
  50. // and only if it is successful.  If it is successful, "plot" may be called
  51. // to actually generate the plot.  Its parameters are as follow:
  52. //
  53. //           f -- z=f(x,y), the function to be plotted.  Before the plot is
  54. //      tilted or rotated, the z-axis runs from the bottom to the top of the
  55. //      display, the y-axis runs from the left to the right of the display,
  56. //      and the x-axis runs out of the display.
  57. //
  58. //           x_min -- the minimum value of x to be plotted.
  59. //
  60. //           x_max -- the maximum value of x to be plotted.
  61. //
  62. //           y_min -- the minimum value of y to be plotted.
  63. //
  64. //           y_max -- the maximum value of y to be plotted.
  65. //
  66. //           external_to_plot -- a function that returns TRUE if and only if a
  67. //      point should be omitted from the plot.
  68. //
  69. //           red -- a function that returns TRUE if and only if a point should
  70. //      be flagged for highlighting.  A point should be so flagged only if it
  71. //      can be seen in the final plot.
  72. //
  73. //           x_division_count -- the number of x divisions to be used in
  74. //      constructing the plot.  At least two must be specified.
  75. //
  76. //           y_division_count -- the number of y divisions to be used in
  77. //      constructing the plot.  At least two must be specified.
  78. //
  79. //           rotation_in_degrees -- rotation (degrees) about an axis parallel to
  80. //      the z-axis and through the center of the surface.
  81. //
  82. //           tilt_in_degrees -- tilt (degrees) about an axis through the center
  83. //      of the surface and parallel to a line from the lower left hand corner of
  84. //      the display to the lower right hand corner of the display.  The plot is
  85. //      tilted after it is rotated.
  86. //
  87. //           (light_x,light_y,light_z) -- a vector pointing to the light source
  88. //      (at infinity).  The light source remains fixed while the plot is rotated
  89. //      or tilted.
  90.     {
  91.       int       display_ready;
  92.       prime_rec initialized_prime_rec;
  93.       int       result;
  94.       int       x_byte_num;
  95.  
  96.       plot_prepared=FALSE;
  97.       if (prime_array_allocated)
  98.         {
  99.           delete prime_array;
  100.           prime_array_allocated=FALSE;
  101.         }
  102.       num_x_divisions=x_division_count;
  103.       num_y_divisions=y_division_count;
  104.       num_primes=(long) num_x_divisions;
  105.       num_primes*=((long) num_y_divisions);
  106.         // Number of quadrilaterals composing the plot.
  107.       initialized_prime_rec.base_z=(unsigned char) '\0';
  108.       initialized_prime_rec.color=(unsigned char) '\0';
  109.       initialized_prime_rec.x=(float) 0.0;
  110.       initialized_prime_rec.y=(float) 0.0;
  111.       initialized_prime_rec.z=(float) 0.0;
  112.       initialized_prime_rec.x_division_index=0;
  113.       initialized_prime_rec.y_division_index=0;
  114.       prime_array=new varray<prime_rec>(initialized_prime_rec,num_primes,5);
  115.         // Virtual array of quadrilaterals composing the plot.
  116.       if (prime_array_allocated=(prime_array->allocated()))
  117.         {
  118.           titillator_ptr=new titillator;
  119.           rotation=rotation_in_degrees;
  120.           tilt=tilt_in_degrees;
  121.           light.x=light_x;
  122.           light.y=light_y;
  123.           light.z=light_z;
  124.           evaluate_and_transform(f,x_min,x_max,y_min,y_max,
  125.            num_x_divisions,num_y_divisions,rotation,tilt,
  126.            external_to_plot,red);
  127.             // Compute the vertices, etc. of the quadrilaterals composing the
  128.             // plot.
  129.           shade();
  130.             // Compute the shade of gray for each quadrilateral composing the
  131.             // plot.
  132.           adjust_perspective();
  133.             // Force parallel lines running away from the viewer to converge at
  134.             // the horizon.
  135.           sort_back_to_front();
  136.           plot_prepared=TRUE;
  137.           delete titillator_ptr;
  138.           display_ready=display_initialized();
  139.             // Do whatever is necessary to prepare the display.
  140.         }
  141.       return (prime_array_allocated & display_ready);
  142.     }
  143.  
  144. void plot3d::evaluate_and_transform(
  145.   double (*f)(double,double),
  146.   double x_min,
  147.   double x_max,
  148.   double y_min,
  149.   double y_max,
  150.   int    num_x_divisions,
  151.   int    num_y_divisions,
  152.   double rotation,
  153.   double tilt,
  154.   int    (*external_to_plot)(double,double),
  155.   int    (*red)(double,double))
  156. // Compute the vertices, etc. for each quadrilateral composing the plot.
  157.     {
  158.       double    cos_rotation;
  159.       double    cos_tilt;
  160.       double    degrees_per_radian;
  161.       double    magnitude;
  162.       prime_rec *prime;
  163.       long      prime_num;
  164.       double    radians;
  165.       double    sin_rotation;
  166.       double    sin_tilt;
  167.       double    tem_x;
  168.       double    tem_y;
  169.       double    tem_z;
  170.       double    x;
  171.       double    x_delta;
  172.       int       x_division_num;
  173.       double    x_rotated;
  174.       double    y;
  175.       double    y_delta;
  176.       int       y_division_num;
  177.       double    z;
  178.  
  179.       degrees_per_radian=45.0/atan(1.0);
  180.       radians=tilt/degrees_per_radian;
  181.       cos_tilt=cos(radians);
  182.       sin_tilt=sin(radians);
  183.       radians=rotation/degrees_per_radian;
  184.       cos_rotation=cos(radians);
  185.       sin_rotation=sin(radians);
  186.       z=f(x_min,y_min);
  187.       x_rotated=x_min*cos_rotation+y_min*sin_rotation;
  188.       y_prime_min=-x_min*sin_rotation+y_min*cos_rotation;
  189.       z_prime_min=-x_rotated*sin_tilt+z*cos_tilt;
  190.       y_prime_max=y_prime_min;
  191.       z_prime_max=z_prime_min;
  192.       x_prime_max=x_rotated*cos_tilt+z*sin_tilt;
  193.       x_delta=(double) (num_x_divisions-1);
  194.       x_delta=(x_max-x_min)/x_delta;
  195.       y_delta=(double) (num_y_divisions-1);
  196.       y_delta=(y_max-y_min)/y_delta;
  197.       x=x_min;
  198.       prime_num=(long) 0;
  199.       for (x_division_num=0; x_division_num < num_x_divisions; x_division_num++)
  200.         {
  201.           titillator_ptr->titillate();
  202.           y=y_min;
  203.           for (y_division_num=0; y_division_num < num_y_divisions;
  204.            y_division_num++)
  205.             {
  206.               z=f(x,y);
  207.               prime=prime_array->vm_ptr(prime_num);
  208.               if (external_to_plot(x,y))
  209.                 prime->base_z=(unsigned char) 3;
  210.               else
  211.                 if (red(x,y))
  212.                   prime->base_z=(unsigned char) 2;
  213.                 else
  214.                   prime->base_z=(unsigned char) 1;
  215.               prime->x_division_index=x_division_num;
  216.               prime->y_division_index=y_division_num;
  217.               x_rotated=x*cos_rotation+y*sin_rotation;
  218.               tem_y=(-x*sin_rotation+y*cos_rotation);
  219.               prime->y=(float) tem_y;
  220.               tem_x=(x_rotated*cos_tilt+z*sin_tilt);
  221.               prime->x=(float) tem_x;
  222.               tem_z=(-x_rotated*sin_tilt+z*cos_tilt);
  223.               prime->z=(float) tem_z;
  224.               if (tem_x > x_prime_max)
  225.                 x_prime_max=tem_x;
  226.               if (tem_y < y_prime_min)
  227.                 y_prime_min=tem_y;
  228.               if (tem_y > y_prime_max)
  229.                 y_prime_max=tem_y;
  230.               if (tem_z < z_prime_min)
  231.                 z_prime_min=tem_z;
  232.               if (tem_z > z_prime_max)
  233.                 z_prime_max=tem_z;
  234.               y+=y_delta;
  235.               prime_num++;
  236.             }
  237.           x+=x_delta;
  238.         }
  239.       magnitude=light.x*light.x+light.y*light.y+light.z*light.z;
  240.       magnitude=sqrt(magnitude);
  241.       light.x=light.x/magnitude;
  242.       light.y=light.y/magnitude;
  243.       light.z=light.z/magnitude;
  244.       return;
  245.     }
  246.  
  247. void plot3d::shade()
  248. // Compute the shade of gray for each quadrilateral composing the plot.
  249.     {
  250.       double     magnitude;
  251.       vertex_rec normal;
  252.       prime_rec  *prime;
  253.       long       prime_num;
  254.       vertex_rec vertex [4];
  255.       int        x_division_num;
  256.       int        y_division_num;
  257.  
  258.       color_min=(unsigned char) (NUM_COLORS-1);
  259.       color_max=(unsigned char) '\0';
  260.       prime_num=num_primes;
  261.       for (x_division_num=num_x_divisions-1; x_division_num >= 0;
  262.        --x_division_num)
  263.         {
  264.           titillator_ptr->titillate();
  265.           for (y_division_num=num_y_divisions-1; y_division_num >= 0;
  266.            --y_division_num)
  267.             {
  268.               prime_num--;
  269.               prime=prime_array->vm_ptr(prime_num);
  270.               vertex[0].x=(double) (prime->x);
  271.               vertex[0].y=(double) (prime->y);
  272.               vertex[0].z=(double) (prime->z);
  273.               if (x_division_num < (num_x_divisions-1))
  274.                 if (y_division_num < (num_y_divisions-1))
  275.                   {
  276.                     prime_num+=((long) num_y_divisions);
  277.                     prime=prime_array->vm_ptr(prime_num);
  278.                     vertex[1].x=(double) (prime->x);
  279.                     vertex[1].y=(double) (prime->y);
  280.                     vertex[1].z=(double) (prime->z);
  281.                     prime_num++;
  282.                     prime=prime_array->vm_ptr(prime_num);
  283.                     vertex[2].x=(double) (prime->x);
  284.                     vertex[2].y=(double) (prime->y);
  285.                     vertex[2].z=(double) (prime->z);
  286.                     prime_num-=((long) num_y_divisions);
  287.                     prime=prime_array->vm_ptr(prime_num);
  288.                     vertex[3].x=(double) (prime->x);
  289.                     vertex[3].y=(double) (prime->y);
  290.                     vertex[3].z=(double) (prime->z);
  291.                     prime_num--;
  292.                   }
  293.                 else
  294.                   {
  295.                     prime_num--;
  296.                     prime=prime_array->vm_ptr(prime_num);
  297.                     vertex[1].x=(double) (prime->x);
  298.                     vertex[1].y=(double) (prime->y);
  299.                     vertex[1].z=(double) (prime->z);
  300.                     prime_num+=((long) num_y_divisions);
  301.                     prime=prime_array->vm_ptr(prime_num);
  302.                     vertex[2].x=(double) (prime->x);
  303.                     vertex[2].y=(double) (prime->y);
  304.                     vertex[2].z=(double) (prime->z);
  305.                     prime_num++;
  306.                     prime=prime_array->vm_ptr(prime_num);
  307.                     vertex[3].x=(double) (prime->x);
  308.                     vertex[3].y=(double) (prime->y);
  309.                     vertex[3].z=(double) (prime->z);
  310.                     prime_num-=((long) num_y_divisions);
  311.                   }
  312.               else
  313.                 if (y_division_num < (num_y_divisions-1))
  314.                   {
  315.                     prime_num++;
  316.                     prime=prime_array->vm_ptr(prime_num);
  317.                     vertex[1].x=(double) (prime->x);
  318.                     vertex[1].y=(double) (prime->y);
  319.                     vertex[1].z=(double) (prime->z);
  320.                     prime_num-=((long) num_y_divisions);
  321.                     prime=prime_array->vm_ptr(prime_num);
  322.                     vertex[2].x=(double) (prime->x);
  323.                     vertex[2].y=(double) (prime->y);
  324.                     vertex[2].z=(double) (prime->z);
  325.                     prime_num--;
  326.                     prime=prime_array->vm_ptr(prime_num);
  327.                     vertex[3].x=(double) (prime->x);
  328.                     vertex[3].y=(double) (prime->y);
  329.                     vertex[3].z=(double) (prime->z);
  330.                     prime_num+=((long) num_y_divisions);
  331.                   }
  332.                 else
  333.                   {
  334.                     prime_num-=((long) num_y_divisions);
  335.                     prime=prime_array->vm_ptr(prime_num);
  336.                     vertex[1].x=(double) (prime->x);
  337.                     vertex[1].y=(double) (prime->y);
  338.                     vertex[1].z=(double) (prime->z);
  339.                     prime_num--;
  340.                     prime=prime_array->vm_ptr(prime_num);
  341.                     vertex[2].x=(double) (prime->x);
  342.                     vertex[2].y=(double) (prime->y);
  343.                     vertex[2].z=(double) (prime->z);
  344.                     prime_num+=((long) num_y_divisions);
  345.                     prime=prime_array->vm_ptr(prime_num);
  346.                     vertex[3].x=(double) (prime->x);
  347.                     vertex[3].y=(double) (prime->y);
  348.                     vertex[3].z=(double) (prime->z);
  349.                     prime_num++;
  350.                   }
  351.               // Compute the normal to a quadrilateral by averaging the
  352.               // normals to each vertex of the quadrilateral.
  353.               normal.x
  354.                =(vertex[1].y-vertex[0].y)*(vertex[3].z-vertex[0].z)
  355.                -(vertex[3].y-vertex[0].y)*(vertex[1].z-vertex[0].z)
  356.                +(vertex[2].y-vertex[1].y)*(vertex[0].z-vertex[1].z)
  357.                -(vertex[0].y-vertex[1].y)*(vertex[2].z-vertex[1].z)
  358.                +(vertex[3].y-vertex[2].y)*(vertex[1].z-vertex[2].z)
  359.                -(vertex[1].y-vertex[2].y)*(vertex[3].z-vertex[2].z)
  360.                +(vertex[0].y-vertex[3].y)*(vertex[2].z-vertex[3].z)
  361.                -(vertex[2].y-vertex[3].y)*(vertex[0].z-vertex[3].z);
  362.               normal.y
  363.                =(vertex[3].x-vertex[0].x)*(vertex[1].z-vertex[0].z)
  364.                -(vertex[1].x-vertex[0].x)*(vertex[3].z-vertex[0].z)
  365.                +(vertex[0].x-vertex[1].x)*(vertex[2].z-vertex[1].z)
  366.                -(vertex[2].x-vertex[1].x)*(vertex[0].z-vertex[1].z)
  367.                +(vertex[1].x-vertex[2].x)*(vertex[3].z-vertex[2].z)
  368.                -(vertex[3].x-vertex[2].x)*(vertex[1].z-vertex[2].z)
  369.                +(vertex[2].x-vertex[3].x)*(vertex[0].z-vertex[3].z)
  370.                -(vertex[0].x-vertex[3].x)*(vertex[2].z-vertex[3].z);
  371.               normal.z
  372.                =(vertex[1].x-vertex[0].x)*(vertex[3].y-vertex[0].y)
  373.                -(vertex[3].x-vertex[0].x)*(vertex[1].y-vertex[0].y)
  374.                +(vertex[2].x-vertex[1].x)*(vertex[0].y-vertex[1].y)
  375.                -(vertex[0].x-vertex[1].x)*(vertex[2].y-vertex[1].y)
  376.                +(vertex[3].x-vertex[2].x)*(vertex[1].y-vertex[2].y)
  377.                -(vertex[1].x-vertex[2].x)*(vertex[3].y-vertex[2].y)
  378.                +(vertex[0].x-vertex[3].x)*(vertex[2].y-vertex[3].y)
  379.                -(vertex[2].x-vertex[3].x)*(vertex[0].y-vertex[3].y);
  380.               prime=prime_array->vm_ptr(prime_num);
  381.               magnitude
  382.                =sqrt(normal.x*normal.x+normal.y*normal.y+normal.z*normal.z);
  383.               if (magnitude == 0.0)
  384.                 {
  385.                   color_min=(unsigned char) '\0';
  386.                   prime->color=(unsigned char) '\0';
  387.                 }
  388.               else
  389.                 {
  390.                   prime->color
  391.                    =(unsigned char) int((double(NUM_COLORS)/2.0)
  392.                    *(1.0+(light.x*normal.x+light.y*normal.y
  393.                    +light.z*normal.z)/magnitude)); // shadows not absolute
  394.                   if (prime->color
  395.                    >= (unsigned char) (NUM_COLORS))
  396.                     prime->color=(unsigned char) (NUM_COLORS-1);
  397.                   if ((prime->color) < color_min)
  398.                     color_min=(prime->color);
  399.                   if ((prime->color) > color_max)
  400.                     color_max=(prime->color);
  401.                 }
  402.             }
  403.         }
  404.       return;
  405.     }
  406.  
  407. void plot3d::adjust_perspective()
  408. // Make parallel lines running away from the viewer appear to converge at the
  409. // horizon.  
  410.     {
  411.       prime_rec  *prime;
  412.       long       prime_num;
  413.       double     tem;
  414.       vertex_rec vertex [4];
  415.       int        x_division_num;
  416.       double     x_eye;
  417.       double     y_center;
  418.       int        y_division_num;
  419.       double     z_center;
  420.  
  421.       if ((y_prime_max-y_prime_min) > (z_prime_max-z_prime_min))
  422.         x_eye=1.1*(y_prime_max-y_prime_min)+x_prime_max;
  423.       else
  424.         x_eye=1.1*(z_prime_max-z_prime_min)+x_prime_max;
  425.       if (((y_prime_max-y_prime_min) > (z_prime_max-z_prime_min))
  426.       ||  (z_prime_max != z_prime_min))
  427.         {
  428.           y_center=(y_prime_max+y_prime_min)/2.0;
  429.           z_center=(z_prime_max+z_prime_min)/2.0;
  430.           prime_num=(long) 0;
  431.           for (x_division_num=0; x_division_num < num_x_divisions;
  432.            x_division_num++)
  433.             {
  434.               titillator_ptr->titillate();
  435.               for (y_division_num=0; y_division_num < num_y_divisions;
  436.                y_division_num++)
  437.                 {
  438.                   prime=prime_array->vm_ptr(prime_num);
  439.                   vertex[0].x=(double) (prime->x);
  440.                   vertex[0].y=(double) (prime->y);
  441.                   vertex[0].z=(double) (prime->z);
  442.                   if (x_division_num < (num_x_divisions-1))
  443.                     if (y_division_num < (num_y_divisions-1))
  444.                       {
  445.                         prime_num+=((long) num_y_divisions);
  446.                         prime=prime_array->vm_ptr(prime_num);
  447.                         vertex[1].x=(double) (prime->x);
  448.                         vertex[1].y=(double) (prime->y);
  449.                         vertex[1].z=(double) (prime->z);
  450.                         prime_num++;
  451.                         prime=prime_array->vm_ptr(prime_num);
  452.                         vertex[2].x=(double) (prime->x);
  453.                         vertex[2].y=(double) (prime->y);
  454.                         vertex[2].z=(double) (prime->z);
  455.                         prime_num-=((long) num_y_divisions);
  456.                         prime=prime_array->vm_ptr(prime_num);
  457.                         vertex[3].x=(double) (prime->x);
  458.                         vertex[3].y=(double) (prime->y);
  459.                         vertex[3].z=(double) (prime->z);
  460.                         prime_num--;
  461.                       }
  462.                     else
  463.                       {
  464.                         prime_num--;
  465.                         prime=prime_array->vm_ptr(prime_num);
  466.                         vertex[1].x=(double) (prime->x);
  467.                         vertex[1].y=(double) (prime->y);
  468.                         vertex[1].z=(double) (prime->z);
  469.                         prime_num+=((long) num_y_divisions);
  470.                         prime=prime_array->vm_ptr(prime_num);
  471.                         vertex[2].x=(double) (prime->x);
  472.                         vertex[2].y=(double) (prime->y);
  473.                         vertex[2].z=(double) (prime->z);
  474.                         prime_num++;
  475.                         prime=prime_array->vm_ptr(prime_num);
  476.                         vertex[3].x=(double) (prime->x);
  477.                         vertex[3].y=(double) (prime->y);
  478.                         vertex[3].z=(double) (prime->z);
  479.                         prime_num-=((long) num_y_divisions);
  480.                       }
  481.                   else
  482.                     if (y_division_num < (num_y_divisions-1))
  483.                       {
  484.                         prime_num++;
  485.                         prime=prime_array->vm_ptr(prime_num);
  486.                         vertex[1].x=(double) (prime->x);
  487.                         vertex[1].y=(double) (prime->y);
  488.                         vertex[1].z=(double) (prime->z);
  489.                         prime_num-=((long) num_y_divisions);
  490.                         prime=prime_array->vm_ptr(prime_num);
  491.                         vertex[2].x=(double) (prime->x);
  492.                         vertex[2].y=(double) (prime->y);
  493.                         vertex[2].z=(double) (prime->z);
  494.                         prime_num--;
  495.                         prime=prime_array->vm_ptr(prime_num);
  496.                         vertex[3].x=(double) (prime->x);
  497.                         vertex[3].y=(double) (prime->y);
  498.                         vertex[3].z=(double) (prime->z);
  499.                         prime_num+=((long) num_y_divisions);
  500.                       }
  501.                     else
  502.                       {
  503.                         prime_num-=((long) num_y_divisions);
  504.                         prime=prime_array->vm_ptr(prime_num);
  505.                         vertex[1].x=(double) (prime->x);
  506.                         vertex[1].y=(double) (prime->y);
  507.                         vertex[1].z=(double) (prime->z);
  508.                         prime_num--;
  509.                         prime=prime_array->vm_ptr(prime_num);
  510.                         vertex[2].x=(double) (prime->x);
  511.                         vertex[2].y=(double) (prime->y);
  512.                         vertex[2].z=(double) (prime->z);
  513.                         prime_num+=((long) num_y_divisions);
  514.                         prime=prime_array->vm_ptr(prime_num);
  515.                         vertex[3].x=(double) (prime->x);
  516.                         vertex[3].y=(double) (prime->y);
  517.                         vertex[3].z=(double) (prime->z);
  518.                         prime_num++;
  519.                       }
  520.                   prime=prime_array->vm_ptr(prime_num);
  521.                   tem=y_center
  522.                    +(vertex[0].y-y_center)*(x_eye-x_prime_max)
  523.                    /(x_eye-vertex[0].x);
  524.                   prime->y=(float) tem;
  525.                   tem=z_center
  526.                    +(vertex[0].z-z_center)*(x_eye-x_prime_max)
  527.                    /(x_eye-vertex[0].x);
  528.                   prime->z=(float) tem;
  529.                   tem=(vertex[0].x+vertex[1].x+vertex[2].x+vertex[3].x)/4.0;
  530.                   prime->x=(float) tem;
  531.                   prime_num++;
  532.                 }
  533.             }
  534.          }
  535.       return;
  536.     }
  537.  
  538. void plot3d::rearrange(
  539.   long lower_bound,
  540.   long upper_bound,
  541.   long *j)
  542.     {
  543.       long      down;
  544.       int       finished;
  545.       prime_rec *prime_1;
  546.       prime_rec *prime_2;
  547.       long      up;
  548.       float     x1;
  549.       int       x_division_index1;
  550.       int       y_division_index1;
  551.  
  552.       prime_1=prime_array->vm_ptr(lower_bound);
  553.       x1=prime_1->x;
  554.       x_division_index1=prime_1->x_division_index;
  555.       y_division_index1=prime_1->y_division_index;
  556.       *j=lower_bound;
  557.       up=upper_bound;
  558.       down=lower_bound;
  559.       do
  560.         {
  561.           finished=FALSE;
  562.           while (! finished)
  563.             if (up <= down)
  564.               finished=TRUE;
  565.             else
  566.               {
  567.                 prime_1=prime_array->vm_ptr(up);
  568.                 if (prime_1->x < x1)
  569.                   finished=TRUE;
  570.                 else
  571.                   up--;
  572.               };
  573.           *j=up;
  574.           if (up != down)
  575.             {
  576.               prime_1=prime_array->vm_ptr(down);
  577.               prime_2=prime_array->vm_ptr(up);
  578.               prime_1->x=prime_2->x;
  579.               prime_1->x_division_index=prime_2->x_division_index;
  580.               prime_1->y_division_index=prime_2->y_division_index;
  581.               finished=FALSE;
  582.               while (! finished)
  583.                 if (down >= up)
  584.                   finished=TRUE;
  585.                 else
  586.                   {
  587.                     prime_1=prime_array->vm_ptr(down);
  588.                     if (prime_1->x > x1)
  589.                       finished=TRUE;
  590.                     else
  591.                       down++;
  592.                   };
  593.               *j=down;
  594.               if (down != up)
  595.                 {
  596.                   prime_1=prime_array->vm_ptr(up);
  597.                   prime_2=prime_array->vm_ptr(down);
  598.                   prime_1->x=prime_2->x;
  599.                   prime_1->x_division_index=prime_2->x_division_index;
  600.                   prime_1->y_division_index=prime_2->y_division_index;
  601.                 }
  602.             }
  603.         }
  604.       while (down != up);
  605.       prime_1=prime_array->vm_ptr(*j);
  606.       prime_1->x=x1;
  607.       prime_1->x_division_index=x_division_index1;
  608.       prime_1->y_division_index=y_division_index1;
  609.       return;
  610.     }
  611.  
  612. void plot3d::quicksort(
  613.   long lower_bound,
  614.   long upper_bound)
  615. //      Recursive quicksort.  According to Knuth, quicksort is suited for
  616. // virtual memory.
  617.     {
  618.       long j;
  619.  
  620.       if (lower_bound < upper_bound)
  621.         {
  622.           titillator_ptr->titillate();
  623.           rearrange(lower_bound,upper_bound,&j);
  624.           quicksort(lower_bound,j-1l);
  625.           quicksort(j+1l,upper_bound);
  626.         }
  627.       return;
  628.     }
  629.  
  630. void plot3d::sort_back_to_front()
  631. //      The painter's algorithm is used; items farther from the viewer are drawn
  632. // earlier.
  633.     {
  634.       quicksort(0l,num_primes-1l);
  635.       return;
  636.     }
  637.  
  638. int plot3d::plot(
  639.   char   *file_name,
  640.   int    show_red,
  641.   int    titillate,
  642.   double bias)
  643. //     This function returns TRUE if and only if it is successful in generating
  644. // the 3D plot.  It calls "aspect_ratio", "pset", and "write_outfile".  Its
  645. // parameters are as follow:
  646. //
  647. //           file_name -- the name of the file to which the plot is to be 
  648. //      written.  For plots on a cathode ray tube, this will usually be "".     
  649. //
  650. //           show_red -- highlight quadrilaterals having each vertex flagged
  651. //      to be highlighted.  Since only those quadrilaterals are redrawn,
  652. //      "plot" should not be called with "show_red" set to TRUE until after 
  653. //      it has been called with "show_red" set to FALSE.
  654. //
  655. //           titillate -- TRUE if the user is to be kept amused while the
  656. //      plot is being generated; FALSE otherwise.  In general, "titillate"
  657. //      should be TRUE for printers and FALSE for cathode ray tubes.
  658. //
  659. //           bias -- a positive number used to adjust the contrast.
  660. //
  661. // "prepare_plot" must be called before "plot", after which "plot" may be called
  662. // as many times as desired.
  663.     {
  664.       double     box_delta_x;
  665.       double     box_delta_y;
  666.       box_rec    box [4];
  667.       int        box_num_1;
  668.       int        box_num_2;
  669.       double     box_x_intercept;
  670.       int        box_x1;
  671.       int        box_x2;
  672.       int        box_y_max;
  673.       int        box_y_min;
  674.       double     box_y_offset;
  675.       int        box_y1;
  676.       int        color_num;
  677.       double     fraction;
  678.       int        intercept_count_mod_2;
  679.       int        line_x1;
  680.       int        line_x2;
  681.       int        outside_maze;
  682.       double     pixels_per_unit;
  683.       prime_rec  *prime;
  684.       long       prime_num;
  685.       int        red_showing;
  686.       int        result;
  687.       int        tint;
  688.       vertex_rec vertex [4];
  689.       int        x_division_num;
  690.       long       x_prime_num;
  691.       int        x_prime_num_mod_50;
  692.       int        y_division_num;
  693.       double     y_offset;
  694.       double     y_out_max;
  695.       double     z_offset;
  696.       double     z_out_max;
  697.  
  698.       if (plot_prepared)
  699.         {
  700.           y_out_max=double(num_x_pixels()-1);
  701.           z_out_max=double(num_y_pixels()-1);
  702.           if (aspect_ratio()*z_out_max*(y_prime_max-y_prime_min)
  703.            > y_out_max*(z_prime_max-z_prime_min))
  704.             {
  705.               pixels_per_unit
  706.                =y_out_max/(aspect_ratio()*(y_prime_max-y_prime_min));
  707.               y_offset=0.0;
  708.               z_offset
  709.                =-(z_out_max-pixels_per_unit*(z_prime_max-z_prime_min))/2.0;
  710.             }
  711.           else
  712.             if (aspect_ratio()*z_out_max*(y_prime_max-y_prime_min)
  713.              < y_out_max*(z_prime_max-z_prime_min))
  714.               {
  715.                 pixels_per_unit=z_out_max/(z_prime_max-z_prime_min);
  716.                 y_offset=(y_out_max
  717.                  -aspect_ratio()*pixels_per_unit*(y_prime_max-y_prime_min))/2.0;
  718.                 z_offset=0.0;
  719.               }
  720.             else
  721.               {
  722.                 pixels_per_unit=1.0;
  723.                 y_offset=y_out_max/2.0;
  724.                 z_offset=-z_out_max/2.0;
  725.               };
  726.           if (titillate)
  727.             {
  728.               x_prime_num_mod_50=0;
  729.               titillator_ptr=new titillator;
  730.             }
  731.           for (x_prime_num=0l; x_prime_num < num_primes; x_prime_num++)
  732.             {
  733.               if (titillate)
  734.                 {
  735.                   x_prime_num_mod_50++;
  736.                   if (x_prime_num_mod_50 >= 50)
  737.                     {
  738.                       x_prime_num_mod_50=0;
  739.                       titillator_ptr->titillate();
  740.                     }
  741.                 }
  742.               prime=prime_array->vm_ptr(x_prime_num);
  743.               x_division_num=prime->x_division_index;
  744.               if (x_division_num < (num_x_divisions-1))
  745.                 {
  746.                   y_division_num=prime->y_division_index;
  747.                   if (y_division_num < (num_y_divisions-1))
  748.                     {
  749.                       prime_num
  750.                        =((long) num_y_divisions)*((long) x_division_num)
  751.                        +((long) y_division_num);
  752.                       prime=prime_array->vm_ptr(prime_num);
  753.                       color_num=(int) (prime->color);
  754.                       if (color_num < NUM_COLORS)
  755.                         {
  756.                           // Adjust contrast.
  757.                           fraction=((double) (color_num-color_min))
  758.                            /((double) (color_max-color_min));
  759.                           if (fraction > 0.0)
  760.                             {
  761.                               fraction=exp(bias*log(fraction));
  762.                               color_num=(int)
  763.                                (((double) (LIGHTEST_GRAY-DARKEST_GRAY))
  764.                                *fraction
  765.                                +((double) DARKEST_GRAY));
  766.                             }
  767.                           else
  768.                             color_num=DARKEST_GRAY;
  769.                         }
  770.                       vertex[0].y=(double) (prime->y);
  771.                       vertex[0].z=(double) (prime->z);
  772.                       red_showing=(prime->base_z == (unsigned char) 2);
  773.                       outside_maze=(prime->base_z == (unsigned char) 3);
  774.                       prime_num+=((long) num_y_divisions);
  775.                       prime=prime_array->vm_ptr(prime_num);
  776.                       vertex[1].y=(double) (prime->y);
  777.                       vertex[1].z=(double) (prime->z);
  778.                       if (red_showing)
  779.                         red_showing=(prime->base_z == (unsigned char) 2);
  780.                       if (outside_maze)
  781.                         outside_maze=(prime->base_z == (unsigned char) 3);
  782.                       prime_num++;
  783.                       prime=prime_array->vm_ptr(prime_num);
  784.                       vertex[2].y=(double) (prime->y);
  785.                       vertex[2].z=(double) (prime->z);
  786.                       if (red_showing)
  787.                         red_showing=(prime->base_z == (unsigned char) 2);
  788.                       if (outside_maze)
  789.                         outside_maze=(prime->base_z == (unsigned char) 3);
  790.                       prime_num-=((long) num_y_divisions);
  791.                       prime=prime_array->vm_ptr(prime_num);
  792.                       vertex[3].y=(double) (prime->y);
  793.                       vertex[3].z=(double) (prime->z);
  794.                       if (red_showing)
  795.                         red_showing=(prime->base_z == (unsigned char) 2);
  796.                       if (outside_maze)
  797.                         outside_maze=(prime->base_z == (unsigned char) 3);
  798.                       if (outside_maze)
  799.                         color_num=NUM_COLORS+1;
  800.                       if ((! show_red) || (red_showing))
  801.                        // Plot each point in a quadrilateral.
  802.                         {
  803.                           if (show_red)
  804.                             color_num=NUM_COLORS;
  805.                           for (box_num_1=0; box_num_1 < 4; box_num_1++)
  806.                             {
  807.                               box[box_num_1].x=(int) (y_offset
  808.                                +pixels_per_unit*aspect_ratio()
  809.                                *(vertex[box_num_1].y-y_prime_min));
  810.                               box[box_num_1].y=(int) (z_offset+z_out_max
  811.                                -pixels_per_unit
  812.                                *(vertex[box_num_1].z-z_prime_min));
  813.                             }
  814.                           box_y_min=box[0].y;
  815.                           box_y_max=box_y_min;
  816.                           for (box_num_1=1; box_num_1 < 4; box_num_1++)
  817.                             {
  818.                               if (box[box_num_1].y < box_y_min)
  819.                                 box_y_min=box[box_num_1].y;
  820.                               if (box[box_num_1].y > box_y_max)
  821.                                 box_y_max=box[box_num_1].y;
  822.                             }
  823.                           for (box_y1=box_y_min; box_y1 <= box_y_max; 
  824.                            ++box_y1)
  825.                             {
  826.                               intercept_count_mod_2=0;
  827.                               box_num_2=1;
  828.                               for (box_num_1=0; box_num_1 < 4; ++box_num_1)
  829.                                 {
  830.                                   if (box[box_num_1].y >= box_y1)
  831.                                     {
  832.                                       if (box_y1 > box[box_num_2].y)
  833.                                         {
  834.                                           box_delta_y=(double)
  835.                                            (box[box_num_2].y
  836.                                            -box[box_num_1].y);
  837.                                           box_delta_x=(double)
  838.                                            (box[box_num_2].x
  839.                                            -box[box_num_1].x);
  840.                                           box_y_offset=(double)
  841.                                            (box_y1-box[box_num_1].y);
  842.                                           box_x_intercept=(double)
  843.                                            (box[box_num_1].x);
  844.                                           box_x1=(int)
  845.                                            ((box_delta_x*box_y_offset)
  846.                                            /box_delta_y+box_x_intercept);
  847.                                           if (intercept_count_mod_2 == 0)
  848.                                             box_x2=box_x1;
  849.                                           else
  850.                                             {
  851.                                               if (box_x1 < box_x2)
  852.                                                 {
  853.                                                   line_x1=box_x1;
  854.                                                   line_x2=box_x2;
  855.                                                 }
  856.                                               else
  857.                                                 {
  858.                                                   line_x1=box_x2;
  859.                                                   line_x2=box_x1;
  860.                                                 }
  861.                                               pset(line_x1,box_y1,
  862.                                                color_num);
  863.                                               while (line_x1 < line_x2)
  864.                                                 {
  865.                                                   line_x1++;
  866.                                                   pset(line_x1,box_y1,
  867.                                                    color_num);
  868.                                                 }
  869.                                             }
  870.                                           intercept_count_mod_2
  871.                                            =1-intercept_count_mod_2;
  872.                                         }
  873.                                     }
  874.                                   else
  875.                                     {
  876.                                       if (box_y1 <= box[box_num_2].y)
  877.                                         {
  878.                                           box_delta_y=(double)
  879.                                            (box[box_num_2].y
  880.                                            -box[box_num_1].y);
  881.                                           box_delta_x=(double)
  882.                                            (box[box_num_2].x
  883.                                            -box[box_num_1].x);
  884.                                           box_y_offset=(double)
  885.                                            (box_y1-box[box_num_1].y);
  886.                                           box_x_intercept=(double)
  887.                                            (box[box_num_1].x);
  888.                                           box_x1=(int)
  889.                                            ((box_delta_x*box_y_offset)
  890.                                            /box_delta_y+box_x_intercept);
  891.                                           if (intercept_count_mod_2 == 0)
  892.                                             box_x2=box_x1;
  893.                                           else
  894.                                             {
  895.                                               if (box_x1 < box_x2)
  896.                                                 {
  897.                                                   line_x1=box_x1;
  898.                                                   line_x2=box_x2;
  899.                                                 }
  900.                                               else
  901.                                                 {
  902.                                                   line_x1=box_x2;
  903.                                                   line_x2=box_x1;
  904.                                                 }
  905.                                               pset(line_x1,box_y1,
  906.                                                color_num);
  907.                                               while (line_x1 < line_x2)
  908.                                                 {
  909.                                                   line_x1++;
  910.                                                   pset(line_x1,box_y1,
  911.                                                    color_num);
  912.                                                 }
  913.                                             }
  914.                                           intercept_count_mod_2
  915.                                            =1-intercept_count_mod_2;
  916.                                         }
  917.                                     }
  918.                                   box_num_2++;
  919.                                   if (box_num_2 >= 4)
  920.                                     box_num_2=0;
  921.                                 }
  922.                             }
  923.                         }
  924.                     }
  925.                 }
  926.             }
  927.           if (titillate)
  928.             delete titillator_ptr;
  929.           result=write_outfile(file_name);
  930.         }
  931.       else
  932.         {
  933.           result=FALSE;
  934.           cerr << "Fatal error:  attempt to plot before prepared." << '\n';
  935.         }
  936.       return result;
  937.     }
  938.